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Thermodynamics of weakly screened (near the one-component-plasma limit) Yukawa fluids in two 
and three dimensions is analyzed in detail. It is shown that the thermal component of the excess in¬ 
ternal energy of these fluids, when expressed in terms of the properly normalized coupling strength, 
exhibits the scaling pertinent to the corresponding one-component-plasma limit (the scalings differ 
considerably between the two- and three-dimensional situations). This provides us with a simple 
and accurate practical tool to estimate thermodynamic properties of weakly screened Yukawa fluids. 
Particular attention is paid to the two-dimensional fluids, for which several important thermody¬ 
namic quantities are calculated to illustrate the application of the approach. 
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I. INTRODUCTION 

Thermodynamics of Yukawa systems (particles inter¬ 
acting via the pairwise Yukawa repulsive potential) is 
of considerable interest, in particular in the context of 
conventional plasmas, dusty (complex) plasmas, and col¬ 
loidal suspensions. For small point-like particles the 
Yukawa potential (also known as screened Coulomb or 
Debye-Hiickel potential) reads 

V(r) = e(A/r) exp(-r/A), (1) 

where e is the energy scale, A is the screening length, 
and r is the distance between a pair of particles. For 
charged particles immersed in a neutralizing medium we 
have eA = Q^, where Q is the particle charge, and the 
screening length A is equal to the corresponding Debye 
radius. Such system are fully characterized by two di¬ 
mensionless parameters. The first is the coupling param¬ 
eter, F = Q^/aT, where a is the characteristic interpar¬ 
ticle separation [in this paper we take a = (47rn/3)“^/^ 
in three dimensions and a = ( 7 rn)~^t '2 jjj dimensions, 
where n is the particle density] and T is the temperature 
(in energy units). The second is the screening parameter, 
defined as k = a/X. 

Statical and dynamical properties of Yukawa systems 
in both two- and three-dimensional conhgurations have 
been extensively investigated using various theoretical 
and simulation techniques. Thermodynamic properties 
have also received considerable attention. The accurate 
results from Monte Carlo (MC) and molecular dynam¬ 
ics (MD) numerical simulations for some thermodynamic 
functions (in particular, for internal energy and com¬ 
pressibility) have been tabulated for a wide (but discrete) 
range of state variables F and k M- Recently, sim¬ 
ple and accurate expressions to evaluate thermodynamic 
properties of three-dimensional (3D) Yukawa fluids Q 
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and solids Q have been proposed, which are very con¬ 
venient for practical applications. These expressions are 
based on the Rosenfeld-Tarazona (RT) scaling Q of the 
excess internal energy in the fluid state (and harmonic 
approximation in the solid state) and deliver impressively 
high accuracy in a wide range of the phase diagram of 
3D Yukawa systems. 

The focus of the present study is on the practical es¬ 
timation of the thermodynamics of Yukawa fluids near 
the one-component-plasma (OCP) limit. The OCP cor¬ 
responds to the limit k —!■ 0 for Yukawa systems. We 
show that for K ;< 1 the OCP scaling of the thermal 
component of the excess internal energy of 3D Yukawa 
fluids is superior to that of the RT scaling. This al¬ 
lows us to somewhat improve the accuracy of our pre¬ 
vious practical ^proach to the thermodynamics of 3D 
Yukawa fluids [3 in the case of weak screening. More 
importantly, we show that the thermal component of the 
excess internal energy exhibits a quasi-universal scaling 
also in two-dimensional systems of soft repulsive parti¬ 
cles (which is, however, quite different compared to the 
RT or OCP scaling in 3D case). Based on the 2D OCP 
limit we then produce a simple analytical expression for 
the excess energy of weakly screened 2D Yukawa systems. 
The accuracy of this expression is verified by comparison 
with the results from MD simulations. We demonstrate 
how this expression can be used to evaluate the specific 
heat, free energy, compressibility factor, and isothermal 
compressibility modulus of 2D Yukawa fluids in practical 
situations. 


II. YUKAWA FLUIDS IN THREE DIMENSIONS 

The reduced excess energy (excess energy per particle 
in units of the system temperature) can be conveniently 
divided into static and thermal components, 

rifl — list “t“ i^th- (2) 

The static contribution corresponds to the value of inter¬ 
nal energy when the particles are frozen in some regular 



2 


structure (e.g. crystalline lattice) and the thermal cor¬ 
rections arise due the deviations from these fixed position 
due to the particle thermal motion. In the strongly cou¬ 
pled classical OCP, the simplest estimate of the static 
component can be obtained using the ion sphere model 
(ISM) [I3,[iil. In this model each particle is restricted to 
the spherical cell of radius a, filled with the neutralizing 
background. Each cell is quasineutral, the cells do not 
overlap, and, therefore, Ust is simply equal to the electro¬ 
static energy of such cell. The latter can be readily cal¬ 
culated, yielding Ugt = m (the energy is negative 

due to particle-background interactions). In fact, this can 
be mathematically proven to be the exact lower bound 
on the excess energy of the OCP The ISM approach 
provides a rather good estimate of the actual OCP excess 
energy at strong coupling, in both fluid and solid phases 
(for example, the Madelung constant of the OCP form¬ 
ing a body-centered-cubic lattice is Mbcc = —0.895929r). 
Another exact lower bound on the OCP excess ener gy i s 
provided by the Debye-Hiickel (DH) construction 
This bound is a reasonable measure of the actual en¬ 
ergy at weak coupling. Recently a physically motivated 
interpolation between the DH and ISM limits has been 
discussed in Ref. M- 

When the value of the static component of the excess 
energy is specified (ISM approach in our case), the ther¬ 
mal component is also well defined. Based on the accu¬ 
rate Monte-Carlo (MC) simulation results from Ref. [l5l| . 
a simple two-term expression has been proposed [T^ 

wth = 0.5944ri/3 - 0.2786. (3) 

This ht, along with the MC numerical data, is plotted in 
Fig- HI The exponent s = 1/3 (or close to that) in Eq. ([3]) 
provides particularly good agreement with the numerical 
data. This fact has been documented in a number of 
previous studies [1, flsl - fisjl . 

For Yukawa systems, the reduced excess energy can 
also be decomposed into the static and thermal compo¬ 
nents [Eq. ([3])], with the static component evaluated us¬ 
ing the generalization of the ISM to these systems (see, 
in particular. Ref. for the details). Note that, in con¬ 
trast to the OCP limit where the neutralizing background 
must be present to keep thermodynamic functions finite 
and well defined, for model Yukawa systems two real¬ 
izations are possible. The first corresponds to the con¬ 
ventional system of charged particles immersed into neu¬ 
tralizing medium, which provides exponential screening 
of the particle-particle electrical interactions and keeps 
the overall charge neutrality (this situation can to some 
extent mimic real colloidal suspensions and complex plas¬ 
mas 0 ). The second is the imaginary one-component 
system of particles interacting via the Yukawa potential 
without any neutralizing medium, which we will refer 
to as the single component Yukawa system (SCYS), and 
which is merely used as a useful model in condensed mat¬ 
ter research. The difference between the excess energies 
of these two systems is related to the interactions be¬ 
tween the particles and the neutralizing medium. The 



FIG. 1: (Color online) Thermal component of the reduced ex¬ 
cess energy, Mth, of the strongly coupled 3D OCP fluid versus 
the coupling parameter P. Crosses correspond to the exact 
data from MC simulations [T^. The (red) solid curve is the 
OCP fit of Eq. ([3]). The (blue) dashed curve corresponds to 
the RT scaling, given by Eq. 0. 


difference can be evaluated as [mill 

3r kT 


Urn = - 


2^2 


( 4 ) 


where the first term represents the energy of the medium 
that, on average, neutralizes the charge of the particles 
(it exactly compensates the energy of the particle-particle 
interactions when the particles are at complete disorder). 
The second term gives the energy of the sheath around 
each particle. The static component of the internal en¬ 
ergy in the ISM approximation is [Til 


, . k{k + l)r kT 3r 
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( 5 ) 


In view of Eq. o, the first term of the above equation 
describes the static energy of the SCYS. The coefficient 
of proportionality, Mg = Ust./T, can be referred to as 
the fluid Madelung constant [^. For SCYS we have 
therefore 


Mfl(K) 


k{k + 1 ) 

(k -I- 1) -I- (/c — 1)62"^ 


( 6 ) 


Note, that the same result [Eq. (jH])] can be obtained using 
the Percus-Yevick (PY) radial distribution function of 
hard spheres in the limit 77 = 1 , where rj is the hard 
sphere packing fraction [ 1 ^ . 

To the best of our knowledge, Rosenfeld and Tara- 
zona were first to discuss a quasi-universal behavior of the 
thermal component of the internal energy (uth = Wex ~ 
MflP) for a wide class of fluids formed by soft repulsive 
particle (in particular, they focused on inversejower-law 
and Yukawa potentials) in three dimensions . The 

accuracy of this scaling for various model systems has 
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been recently investigated in extensive numerical simu¬ 
lations [ 2 ^. It has been observed that the RT scaling 
is particulary useful for liquids with strong correlations 
between equilibrium fluctuations of virial and potential 
energy [ 2 ^, which are referred to as Rosklide-simple or 
just Rosklide systems [ 2 ^. OCR and strongly coupled 
Yukawa systems belong to this class and, thus, RT scal¬ 
ing can be quite useful to describe thermodynamic prop¬ 
erties of these systems. 

Previously, a variant of the RT scaling has been suc¬ 
cessfully used to obtain practical expressions for the in¬ 
ternal energy and pressure of 3D Yukawa fluids with 
K ^ 5 across coupling regimes Q- The expression for 
the thermal component of the excess energy suggested in 
this work is 

uth^3.2(r/r^)2/5-o.i, (7) 

where Tm is the coupling parameter at the fluid-solid 
phase transition. More recently, it has been also pointed 
out that for Yukawa fluids sufficiently close to freezing 
an even simpler scaling, uth — 3.1(r/ri„)^^^, is more 
appropriate [8|. 

The first question to be addressed now is which of the 
scalings, OCP [Eq. ([3])] or RT [Eq. ([7|)] is more appro¬ 
priate in the OCP limit. Figure [1] shows the comparison 
of these scaling with the exact MC simulation data from 
Ref. [ 13 . To plot the RT scaling we have used Pm — 171.8 
for the OCP in three dimensions, as reported in Ref. [ 2 ^ 
(the true Pm for the OCP can be slightly higher, around 
cs 174 [II, [3; we use Pm — 171.8 here for the sake of 
consistency with the Yukawa case, see below). Figure [1] 
demonstrates that the RT scaling describes fairly well the 
numerical data, but the OCP scaling is somewhat more 
accurate, especially in the regime 1 ^ P ^ 100. This is 
absolutely not surprising, since the OCP scaling is noth¬ 
ing but the best simple fit to the numerical data shown 
in Fig. [TJ 

The second relevant question is which of the scalings 
is more relevant for Yukawa systems near the OCP limit 
(i.e., at weak screening). To answer this, we have used 
the numerical data for the excess energy of Yukawa fluids 
tabulated in Ref. subtracted the static contribution 
[Eq. ([5])], and plotted the remaining thermal component 
as a function of P/Pm in Fig. (21 T o do that we used the 
fit for rm(K), proposed in Ref. 

rm(K) 171.8-b 42.46^2-b3.84lK^ (8) 

which is applicable for k < 1.4. The corresponding nu¬ 
merical data for n < 1.0 are shown in Fig. [2] along with 
the OCP and RT scalings. We see that the OCP scaling 
is a better approximation for P/Pm ^ 2 x 10“^. The 
deviations between these two scalings are particularly 
observable for intermediate coupling P ^ 0(10), as in¬ 
dicated in Fig. [2] The observation that the OCP scaling 
provides better description for the thermal component 
of the excess energy of weakly screened (k < 1) Yukawa 
systems allows us to improve the previous practical ap¬ 
proach to the thermodynamics of Yukawa fluids in this 



FIG. 2: (Color online) Thermal component of the reduced ex¬ 
cess energy of weakly screened («;<!) Yukawa fluids in three 
dimensions as a function of the reduced coupling parameter 
r/Fm. Symbols are the MC simulation results from Ref. [^. 
The (red) solid curve shows the OCP scaling based on Eq. Q 
and the (blue) dashed curve corresponds to the RT scaling, 
given by Eq. 0. The vertical lines denote the two groups 
of symbols, corresponding to F = 10 and F = 20, where the 
deviation between these two scalings is most pronounced. 



FIG. 3: (Color online) Reduced excess isochoric heat capac¬ 
ity Cv of the one-component-plasma in 3D versus the coupling 
parameter F. Crosses are the MC results by Brush, Sahlin, 
and Teller (BST) [^. Different crosses for the same F corre¬ 
spond to the systems with different number of particles and 
configurations (for details see Ref. H^). Circles denote the 
MC results by Hansen [^. The (red) solid curve is based 
on the OCP scaling of Eq. 0 and the (blue) dashed curve 
corresponds to the RT scaling of Eq. ©. 


regime. The modifications to the expressions derived in 
Ref. 0 are straightforward and we do not discuss this 
further in detail, except for a relevant comment below. 

It is clear that the thermodynamic functions depend¬ 
ing on the total excess energy (like e.g. pressure and 
compressibility modulus), will be only slightly affected 
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by the differences between the OCP and RT scalings of 
the thermal energy, since the static energy component 
dominates for soft repulsive potentials. However, other 
thermodynamic functions can be more sensitive to these 
differences. For example, the reduced specific heat at 
constant volume is determined by the reduced thermal 
energy alone, 

cv(k, r) = | + Mth-r(5uth/9r). (9) 

In Figure [3] we plot the available numerical data for Cv 
for the OCP, along with the analytical expressions based 
on the OCP and RT scalings. The scattering of the nu¬ 
merical data points does not allow us to clearly discrim¬ 
inate between the OCP and RT scalings. Nevertheless, 
based on the comparison shown in Figs. [T] and [H it is 
not unreasonable to expect that the OCP scaling is more 
appropriate for very soft repulsive potentials, which im¬ 
plies, for instance, instead of ^ y-2/5 

(according to the RT scaling) in this regime. 


III. YUKAWA FLUIDS IN TWO DIMENSIONS 

As in the previous section, we first consider the OCP 
limit. In two-dimensions, two different systems are ac¬ 
tually referred to as the OCP. The first is character¬ 
ized by logarithmic interactions, V(r) = —Q^\n{r/L), 
where L is an arbitrary length scale [28l - [30l |. The loga¬ 
rithmic potential emerges from the solution of the two- 
dimensional Poisson equation around the central point¬ 
like particle [l^. In this case the system is convention¬ 
ally characterized by the coupling parameter F = jT, 
which does not depend on the particle density. In the sec¬ 
ond system, the interaction is of conventional Coulomb- 
type, V{r) = Q^/r, but the interacting particles are con¬ 
fined to the 2D plane. In this second case, the coupling 
parameter is defined similarly to the OCP (or Yukawa) 
system in 3D, F = Q^/aT, with a being the 2D analog 
of the Wigner-Seitz radius (defined in the Introduction). 

The numerical results for the internal energy of both 
OCP systems are available in the literature [291 - 13^ . The 
first interesting point to address, is wether the scaling 
of the thermal component of the excess energy in two- 
dimensions is also quasi-universal for soft repulsive sys¬ 
tems (like the two OCPs considered here), similarly to 
the 3D case. For the static component of the energy in 
the 2D case we use the triangular lattice sums (Madelung 
energies), which are um = — 0.37438r for the logarithmic 
and Um = — 1.106103F for the Coulomb potential, re¬ 
spectively [ 2 ^ 113 . (The ion disc model, an analog of the 
ISM in 3D, can be constructed for logarithmic interac¬ 
tions [ 2 ^ , but we are not aware of any such construction 
for the Coulomb interaction in 2D). Then, subtracting 
the static component from the full excess energy (when 
necessary) we can obtain the thermal energy component. 
The resulting dependence of Uth on the reduced coupling 
parameter F/Fm is shown in Figure SI To produce this 
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FIG. 4: (Color online) Thermal component of the reduced 
excess energy of several soft repulsive systems in two dimen¬ 
sions as a function of the reduced coupling parameter F/Fm- 
Circles correspond to the numerical data for the Coulombic 
interaction in 2D [^, triangles show the MC results for log¬ 
arithmically repelling particles confined to the surface of a 
sphere [ 2 ^, while squares correspond to the MD results for 
logarithmic repulsion on a plane 1301] . The crosses correspond 
to the IPL potential V{r) oc r~" with n = 6 [sj. The (red) 
solid curve shows the OCP scaling for two dimensional sys¬ 
tems, Eq. CQl). The brown curve is the fit to the n = 6 IPL 
data using the same functional form of Eq. (HQl). 

figure we assumed Fm = 140 for the logarithmic inter¬ 
action fl^ [^ . [ 30 II and Fm = 137 for the Coulomb in¬ 
teraction . The closeness of Fm values for these two 
systems is likely a coincidence, since the physical mean¬ 
ing of the coupling parameters is quite different for the 
logarithmic and Coulomb interactions. 

The numerical data for the thermal component of the 
excess energy for 2D OCP systems shown in Fig. 0] tend 
to collapse on a single quasi-universal curve, which is 
analogous to the OCP scaling in three dimensions. Some 
scattering of the data points is present, but no clear sys¬ 
tematic trend is observed, indicating that this can sim¬ 
ply reflect the level of accuracy of the simulation results. 
The dependence of Uth on F/F^ has a logarithmic char¬ 
acter (in contrast to the power-law scaling in 3D) and we 
found convenient to approximate the data by the follow¬ 
ing functional dependence: 

Mth ^iin[i-i- 62 (r/rm)], (10) 

with the coefficients bi = 0.231 and 62 = 391.655. Also 
shown in Fig. |4]are the numerical data for the thermal 
energy of the fluid with inverse power law (IPL) inter¬ 
action, V{r) = e{alrY, with n = 6 [s^. To produce 
this plot we have used the notation F = {e/T){a/a)^, 
and took Fm ~ 105.9 at the fluid-solid phase transition 
(in fact, this is the fluid-hexatic transition) [s^. We see 
that the universal scaling only holds for sufficiently soft 
interactions, which is not the case for the IPL n = 6 case 
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(the thermal energy is systematically lower). However, 
the fit of the form can still be reasonably applied, 
see the dashed curve in Fig. |4] (the fit yields bi = 0.251 
and 62 = 163.290 for the n = 6 IPL). 

Now, by analogy with approach discussed in the previ¬ 
ous section, we apply the observed soft-repulsion scaling 
of Mth to evaluate the excess energy of weakly screened 
Yukawa systems in 2D geometry. The static component 
of the excess energy is associated with the triangular lat¬ 
tice sum (Madelung energy), which for the Yukawa in¬ 
teraction can be approximated by @ 

UM ^ -1.106ir (1 - 0.4555k -b 0.0999^2 - 0.0087 k^) . 

( 11 ) 

Note that this static energy accounts for particle-particle 
correlations and the energy of the neutralizing medium 
(this is why the energy is negative). The energy of the 
neutralizing medium (plasma) compensates the energy 
coming from particle-particle interactions when no cor¬ 
relations are present. It can be therefore calculated us¬ 
ing the virial energy equation with the radial distribution 
function g(r) = 1 , that is 

pOO -p 

/ rV{r)dr = —, (12) 

where V (r) is the Yukawa potential. The term —Upi is 
sometimes referred to as the positive Hartree part, while 
the total energy is then called the correlational part . 
Obviously, the energy of the SCYS in two dimensions can 
be obtained as a sum of the Hartree and the correlational 
parts. (The total energy may also include the contribu¬ 
tion from particle-sheath interactions, similarly to the 3D 
case. However, this term does not affect main thermody¬ 
namic quantities of interest and we drop it from further 
consideration). 

The dependence rm(K) for two dimensional Yukawa 
systems has been approximated in Ref. by the fol¬ 
lowing fit 

p c^') ^_ (13") 

’ 1 - 0.388k2-p0.138k3-0.0138k4' ^ ’ 

This fit underestimates by ~ 4% the transitional coupling 
parameter of the OCP, but describes relatively well the 
melting points found from the bond angular correlation 
analysis (see Fig. 6 of Ref. [l^) up to k = 3.0. 

The combination of Eqs. (HU) - dni) allows us to cal¬ 
culate the energy of strongly coupled Yukawa fluids in 
2D. We compare our results with the energy data tabu¬ 
lated in Ref. Q, which are in good agreement with the 
subsequent simulations in Ref. [s^. The comparison is 
shown in Figure [5] The agreement is rather good, the 
maximum relative deviation from the numerical results 
does not exceed few percent and is much smaller than 
that for the majority of the data points. 

The data tabulated in Ref. Q do not cover the 
regime of very strongly coupled fluids, near the fluid- 
solid phase transition. We, therefore, performed addi¬ 
tional simulations for 2D Yukawa fluids near freezing. 



FIG. 5: (Color online) Reduced excess energy Wex/F as a func¬ 
tion of the coupling parameter F for weakly screened two 
dimensional Yukawa fluids. The symbols correspond to the 
numerical results from Ref. [^. The solid curves represent 
our approximation. Shown are the results for five screening 
parameters: k = 0.0, 0.5, 1.0, 1.5, and 2.0. 


Simulations have been performed on graphics process¬ 
ing unit (NVIDIA GTX 960) using the hoomd-blue soft¬ 
ware [H, 113 • We used N = 160000 Yukawa particles 
in a rectangular box with the periodic boundary condi¬ 
tions. The cutoff radius for the potential has been chosen 
Lcut — 26.6A. Simulations have been performed in the 
canonical ensemble {NVT) with the Langevin thermo¬ 
stat at a temperature corresponding to the desired target 
coupling parameter F. We first equilibrated the system 
for one million time steps and measured and averaged 
the excess energy for another million of steps. The re¬ 
sults of this simulation along with their comparison with 
the theoretical predictions is summarized in Table HI The 
agreement between the theory and numerical experiment 
is excellent for k < 1.5. In this regime the deviations 
between the theoretical and MD results are all within 
the simulation uncertainty. For k = 2.0, the theory sys¬ 
tematically overestimates the excess energy obtained in 
simulations, but the relative difference is below 1%. This 
qualitative trend is not unexpected. As the potential 
softness decreases, the OCP scaling tends to overestimate 
the actual thermal component of the excess energy. An 
example, corresponding to n = 6 IPL potential, is shown 
in Fig. m 


Summarizing this Section, we have proposed simple 
and reliable practical method to estimate thermodynam¬ 
ics of weakly screened Yukawa fluids in two dimensions. 
The accuracy of this method is very good, at least for 
K 2. In the next sections we provide an exemplary 
calculation of several thermodynamic functions. 
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TABLE I: Excess energy of the two-dimensional single component Yukawa fluid in the vicinity of the fluid-solid transition. 
Three coupling regimes, characterized by T/Tm = 0.8, 0.9, and 0.95 have been investigated [Tm has been evaluated according 
to Eq. (I13ll 1. For each of these regimes four values of the screening parameter k have been chosen (first column). The results 
from numerical MD simulations (MD) are compared with the theoretical prediction (Theory). 


K, 

F/Fn, = 0.80 


r/r,n = 0.90 


r/r,n = 0.95 


MD 

Theory 

MD 

Theory 

MD 

Theory 

0.5 

128.95 ±0.34 

128.93 

144.82 ± 0.35 

144.90 

152.79 ±0.36 

152.89 

1.0 

43.62 ±0.10 

43.59 

48.88 ±0.12 

48.90 

51.55 ±0.13 

51.55 

1.5 

21.42 ±0.05 

21.43 

23.95 ± 0.06 

23.97 

25.24 ± 0.06 

25.24 

2.0 

12.89 ±0.03 

12.99 

14.36 ± 0.03 

14.48 

15.09 ±0.04 

15.22 




FIG. 6: (Color online) Reduced isochoric heat capacity Cv 
of the two-dimensional weakly screened Yukawa fluid as a 
function of the coupling parameter F. The curves denote the 
calculation using the present approximation for four values of 
the screening parameter: k — 0.5, 1.0, 1.5, and 2.0. 


FIG. 7: (Color online) Reduced excess free energy /ex of the 
two-dimensional weakly screened Yukawa fluid as a function 
of the coupling parameter F. The notation is the same as in 

Fig. El 


IV. THERMODYNAMIC FUNCTIONS OF 
TWO-DIMENSIONAL YUKAWA FLUIDS 

We use the conventional thermodynamic identities [s^ 
to calculate several useful quantities for Yukawa fluids in 
two dimensions. In doing so it is convenient to operate 
with the Yukawa system phase state variables k and F. If 
the density and temperature of the neutralizing medium 
are kept fixed, we have F oc {aT)~^ oc and k oc 

a oc This implies 

dr _ F 9F _ 1F dK_ dn _ Ik 

'Wf~~T" W'~ ’ 

The practical model expression for Uey^{K, F) can then be 
employed to calculate various thermodynamic quantities. 

The first example deals with the reduced isochoric heat 
capacity, which is related to the excess energy via Cv = 
2 -I- Uex — F(9uex/9F). Since the static contribution is 
linear in F, Cv is completely determined by the thermal 
component of the excess energy, see Eq. (0. We have 
calculated the dependence of Cv on F for several values 
of K (k < 2.0) and plotted the resulting curves in Fig. |6l 

The second quantity we consider is the reduced excess 
free energy of Yukawa fluids. It is related to the excess 


energy by a standard integration, 

/ex(K,F)= [ dF'Mex(K,F')/F'. (14) 

Jo 

An exact result would be obtained if the dependence 
Mex(K,r) is known in the entire range of coupling, start¬ 
ing from F = 0. This is not the case here, because 
we have proposed an approximation for the sufficiently 
strong coupling regime (say F > 1). The conventional 
procedure can be applied to start the integration from 
F = 1 and add the weak coupling correction /ex(K, 1). 
The latter can be calculated using the analytical results 
available in the literature (see e.g. @). However, this 
contribution is vanishingly small at F 3> 1 and can be 
omitted. Since both u^t and itth are linear in F in the 
weak coupling limit, we will not introduce serious errors 
by starting the integration in Eq. (fHl) from F' = 0. The 
obtained results for the reduced free energy are shown 
in Fig. [T] These compare very favorably with the data 
presented in Fig. 3 of Ref. @ (including the regime F < 1 
shown in the inset of this figure). 

Finally, we calculate the compressibility factor (re¬ 
duced pressure). 


Z (k, F) — 1 -|- 


Tdu 


2 9F 


K dfe^ 
2 dK ’ 


(15) 
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FIG. 8: (Color online) Compressibility factor Z (a) and 
the isothermal compressibility modulus [i (b) of the two- 
dimensional weakly screened Yukawa fluid as functions of the 
coupling parameter F. The notation is the same as in Fig. (6] 

and the isothermal compressibility modulus, 

V) = Z + {T/2){dZ/dT) - (k/2)(5Z/5k). (16) 

These quantities can be of considerable interest in the 
context of wave phenomena in Yukawa fluids [s^. The 
calculated Z and /i, as functions of F, are plotted in Fig. [8] 
for several values of the screening parameter k. Our 
present calculation for the compressibility factor com¬ 
pares very favorably with the MD simulation results pre¬ 


sented in Fig. 4 of Ref. @ (again, including the regime 
F ;< 1 shown in the inset of this figure). Thus, simple 
and accurate practical tool to estimate thermodynamic 
properties of 2D Yukawa fluids is now available. 


V. CONCLUSION 

The Rosenfeld-Tarazona (RT) scaling i 113 of the 
thermal component of the excess energy of simple fluids 
with soft repulsive interactions is a useful tool to esti¬ 
mate their thermodynamic properties. Previously, the 
approach based on the RT scaling arguments was ap¬ 
plied to three-dimensional Yukawa fluids and very good 
accuracy was documented [3, 0 • this paper we have 
shown that there is still room for some slight improve¬ 
ments in the regime of weak screening. In this regime, 
the one-component-plasma-like scaling of the thermal en¬ 
ergy is superior to the RT scaling and thus even higher 
accuracy in estimating thermodynamic properties can be 
reached. Moreover, we have documented the existence of 
similar scaling (quasi-universal behavior) of the thermal 
energy for two-dimensional fluids with soft repulsive in¬ 
teractions (the functional form of the scaling varies sig¬ 
nificantly between the two- and three-dimensional cases). 
Based on that the total excess energy of weakly screened 
two-dimensional Yukawa fluids has been evaluated and 
good agreement with both the existing and new results 
from direct molecular dynamics simulations has been ob¬ 
served. To further demonstrate the application of this 
approach, various important thermodynamic quantities 
of the two-dimensional Yukawa fluids have been evalu¬ 
ated in a wide range of coupling. 
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